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ΟΙ АД 


The complex geometry of many parts of a ceramic gas tur- 
bine engine and the need to know the stresses in these parts 
to a high degree of accuracy indicate that the finite element 
method should be employed. One analysis code used a finite 
element program based on an isoparametric 12-node brick. 

This allowed quadratic variation in one coordinate direction, 
but only linear variation in the other two coordinate direc- 
tions. This introduced an artificial stiffness into the 
structure. This study was made to determine the effects of 
uSing a totally quadratic isoparametric brick, vice one 
which is quadratic in only one direction. 

A distributed load preprocessor was also developed and 
tested. | 


FORTRAN IV was used throughout. 
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Е И ООО ОЕ 


ΠΥ CERAMICS? 

There are three major advantages in using ceramics in 
the hot section of a gas turbine: improved engine perform- 
ance, higher durability, and less dependence on expensive 
imported materials. 

Simple Brayton cycle analysis shows that a higher tur- 
bine inlet temperature (TIT) results in a higher thermal 
efficiency, and lower specific fuel consumption. Blade 
Mie terval considerations usually limit TIT. Intricate blade 
cooling schemes, such as those used in the 1M2500 and the 
FT9, use impingement and convective heát transfer and exter- 
nal film cooling, to lower the temperature of the turbine 
blade. This degrades engine performance and introduces 
manufacturing complexities. Additionally, the cooling 
holes are susceptible to plugging, which may lead to even- 
tual failure. 

Ceramics, however, can withstand higher turbine inlet 
temperatures without blade cooling. Ceramic components 
also have significantly higher resistance to hot corrosion 
than metal "superalloy" engines. The corrosion properties 
of ceramics make them better candidates to withstand the 
higher fuel impurities expected in the future. 

At present, metal engines require expensive, imported 


metals such as chromium, cobalt, columbium, and nickel. 
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These metals are in limited supply and come from areas of 
potential political instability. Ceramics rely on inexpen- 


πο αυ materials which are found in the United States. 


S DESIGN CONSIDERATIONS 

Ceramics are brittle materials. This does not mean 
that they are necessarily fragile. A brittle material will 
fracture with little or no plastic deformation. A fragile 
material will fracture at low stress with little or no 
plastic deformation. It has been suggested that the word 
"inductile" be used instead of brittle [1]. The inductile 
nature of ceramics is the crux of the design problem with 
this material. Local stress concentrators, such as the tip 
of a crack, will not yield in plastic deformation. If the 
stress concentration is not relieved, material fracture m 
pesult. 

Another major problem is the variability of data on the 
material properties of ceramics. 

These two problems have a major impact on the design. 
Design procedures that were valid for most engineering 


metals may fail when applied to ceramics. 


EN PROJECT BACKGROUND 

Some of the calculations performed later in this work 
are concerned with the design of a ceramic engine for a 
Defense Advanced Research Projects Agency (DARPA) project. 
The prime contractor, the AiResearch Division of the 


Garrett Corporation, used their T76 model metal engine as 
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the basis for the ceramic design. They indicated that by 

replacing certain metal hot section parts, they could im- 

prove engine output power by 40%, and reduce specific fuel 
consumption by 10%, as shown in Fig. 1. An engine cross- 

section (Fig. 2) shows the ceramic hot section. 

As mentioned earlier, the brittle nature of ceramics 
drives the design. In the case of the Garrett engine, a 
compliant layer was added in the contact region of the 
dovetail as an attempt to reduce the contact stresses. 
Blade tip clearance was increased to reduce the possibility 
of impact of the blade on the shroud. In addition, a prob- 
abilistic approach, based on the Weibull distribution, 
was used as part of a design methodology that was aimed at 
predicting probability of success. These are but a few of 


the examples of how ceramics affected the design. 


ШИ REASON FOR THESIS 

The complex geometry of many of the ceramic engine's 
parts, and the need to know the stresses to a high degree 
of accuracy, indicated that the finite element method 
Should be employed. One analysis code used a finite element 
program based on an isoparametric 12-node brick. This 
allowed quadratic variation in one coordinate direction, 
but only linear variation in the other two coordinate direc- 
tions. This introduced an artificial stiffness into the 
structure. This study was made to determine the effects of 
using a totally quadratic isoparametric brick, vice one 


Which is quadratic in only one direction. 
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TO HEOL Y 


Recall from strength of materials: 


M = Emp 


where M = moment, E = Young's modulus, I = moment of iner- 
tia, and p = radius of curvature. 

Consider a beam made of an isotropic material, with a 
constant cross-section. If the beam is in pure bending, 
the moment is constant across the span. If Young's modulus 
and moment of inertia are also constant, the radius of cur- 
vature is a constant. A curve which has a, constant radius 
Of curvature is an arc of a circle. 

Figure 3(a) shows a 4-node plane element. There can 
only be linear variation on each of the sides. If the ele- 
ment is deformed in pure bending, each of the sides remains 
straight. The top and bottom surfaces do not reflect any 
of the curvature expected. This is a very stiff approxima- 
Оп. 

A 6-node plane element, shown in Figs. 3(b) and 3(c), 
yields a better approximation to the bending phenomenon, 
depending on the orientation of the one quadratic side. 
Orientation 1 does not reflect any of the bending curvature. 
Orientation 2, however, does. The finite element method 
would represent the deflected curve, which is defined by 3 


nodes, by passing a parabola through the node points. This 
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can give a very close approximation to the circular are 
predicted by strength of materials. It should be noted, 
however, that orientation 2 is constrained by the same 
"planes remain plane" restriction applied to slender beam 
theory. Thus, depending on the geometry and loading, a 6- 
node plane element with orientation 2 may or may not accu- 
rately reflect bending for a deep beam. 

A totally quadratic element, shown in Fig. 3(d), is 
not constrained by orientation, geometry or loading. An 8- 
node plane element can give a better representation of bend- 
ing. 

For a given material, a stiffer element, such as the 4- 
or 6-node plane element, will deform less than it should. 
Since strain is based on displacement, and stress is based 
Ae train, the stiffer structure will predict a lower stress 
than is actually present. This is a potentially dangerous 
design situation. 

The development presented here is based on a two-dimen- 
sional plane element. Similar reasoning can be applied to 
three-dimensional bricks. A 12-node brick, for example, 
has only one quadratic side. Depending on the loading, 
geometry, and orientation of the quadratic side, a 12-node 
brick could approach the correct solution. It will always, 
however, be stiffer than the actual structure. If it is 
oriented incorrectly with respect to the applied moment, it 
could result in seriously incorrect answers for a coarse 


mesh system. In addition, since two of the coordinate 


LO 





directions are stiffer than the third quadratic coordinate 
direction, an artificial anisotropy is introduced into the 
тореш. 

A fully quadratic brick does not have the problems men- 
tioned above. Although the discretization process intro- 
duces some artificial stiffness, it is very slight. A 
fully cubic brick would even be a better representation of 


mne actual structure. 


17 





ME NDINE SEE CEMENTOS LED 


To evaluate the accuracy of a finite element computer 
program, it is often useful to compare its results to experi- 
mental results for a problem that can be verified experimen- 
tally. A 4-point flexure test has been used very often for 
this purpose [3]. This chapter explains how the flexure 
test can be modeled using 12-node and 20-node bricks. The 


K rS are also discussed. 


СЕИ РООЕ POINT FLEXURE TESTS 

Figure 4 shows a schematic of a 4-point flexure rig. 
It is used to provide test results from a specimen in pure 
bending. Figure (5a) shows a free body diagram for the 
specimen. Figures (5b) and (5c) show the shear and moment 
diagrams. The moment diagram shows that the portion of the 
span between the applied forces on the top surface has no 


shear, and maximum constant moment, i.e., pure bending. 


B. MODELING TECHNIQUES 

ο... πο νο балон x 0.250 inch x 3.0 inch bend 
specimen was used. The specimen was modeled using a vari- 
ably-noded three-dimension brick from the ADINA [4] and 
SAP IV [5] finite element programs. Maximum use was made 
of the symmetry in the problem. Only 25% of the specimen 
needed to be modeled, as shown in Figs. 6a,b,c. Symmetry 


in the spanwise and long tranverse directions was used. 
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A cross-sectional view of the deformed bar, Fig. 7, shows 
the symmetry in the long transverse direction. The mid-span 
cross-section was given a guided end condition (no displace- 
ment in the x-direction). The midplane in the long trans- 
verse direction was also given a guided end condition (no 
displacement in the y-direction). 

Since the model was very regular in shape, automatic 
nodal point coordinate and connectivity generation were ex- 
tensively used. Both were generated in the spanwise direc- 
tion to take maximum advantage of the automatic generation 


features of the codes employed. 


ιν. LEXAURE CASE STUDY 

А 12-поде brick has one quadratic direction. The other 
two directions are linear. When used to model the bending 
specimen, three distinct orientations of the brick can occur. 
Figure 8 shows each of these orientations with respect to 
the bending moment. An analysis was performed using the 
SAP IV program for each of the three orientations, called 
Case 1, Case 2, and Case 3. The applied forces were modeled 
as distributed loads. 

1. Analysis of Results 

The three-dimensional finite elements used to solve 

the problem would give a solution converging to an exact 
three-dimensional elasticity theory solution. However, such 
a solution in closed form does not exist. Therefore, the 
results were compared to the usual strength of materials 


beam theory. Such a theory is not perfect and circumspection 
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had to be used in such a comparison. For instance, beam 
theory does not predict any anticlastic curvature [18], but 
the finite element solution shows some. However, due to 
the boundary conditions and symmetry, beam theory should 
apply to the x-z plane of symmetry with a very good accura- 
су. The displacements of the neutral axis of this plane 
from x = 0 to x = L were plotted versus the span length. 
Displacements were normalized by dividing them by the mid- 
span (x = 0) displacement predicted by beam theory (53). 
This facilitates comparisons in terms of the percentage of 
the maximum deflection of beam theory. 

Two loadings and materials were used. The first 
Specimen was made of steel, and had a 300 1b load. The 
second specimen was hot pressed silicon nitride (HPSN), 
and had a 180 1b load. The mechanical properties of HPSN 
are given in Ref. [6]. 


The deflection curves were calculated using [7]: 


2 
ο ES _ ο. 
Α СЕТ СОЕ 
_ w(L-a) 2 2 
where Sa MEET ze ы 2aL - a) 
ПИ пен = reaction end moment 


w = applied force 
E = Young's modulus 


I enementr er cmertıa 


Figure 9 shows that L is half the distance between supports 
on the lower side of the specimen, and w is the total load 


ο ήτο, i.e., 300 Ib or 180 ib. The symbol < > denotes a 
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Anseontinmsous funetion, 1.e., for x<a the quantity <x-a> = 0, 
mer x>07°<x-a> = x-a. 

The results of these analyses are shown in Figs. 10 
and 11. The Case 1, 2 and 3 curves correspond to three 
orientations of the 12-node brick. Case 1, the orientation 
with the quadratic side along the axis of bending, is the 
closest of the three to the values predicted by beam theory. 
The maximum deflection of the ceramic specimen is approxi- 
mately 96% of the deflection predicted by beam theory. This 
model is slightly stiffer than beam theory. Cases 2 and 3, 
which have the quadratic side perpendicular to the axis of 
bending, are significantly stiffer models than Case 1. The 
maximum deflection predicted by these cases is only 70% of 
the deflection predicted by beam theory. 

In the displacement-based finite element method, 
stress and displacement are related by [8]: 


{е} πι 


Ciel [O| [pi lu) 


апа (t) 


where (т) = a vector composed of the six independent 
components of the stress tensor. 
{e} = a vector composed of the six independent 
components of the small strain tensor. 
[B] = strain-displacement transformation matrix. 
{u} = displacement matrix. 
[C] = matrix of material properties. 
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This formulation clearly shows that a stiffer struc- 
ture (i.e., less displacement and displacement gradients) 
will produce less stress. 

Case 4, shows results which are slightly more flex- 
ible than beam theory. Since beam theory is inherently 
stiff due to the approximations made during its development 
(i.e., negligible shear stress), the Case 4 model is closer 
to predicting the actual stress distribution than either 
Cases 1, 2, 3 or beam theory. It should also be remembered 
that the actual structure modeled will bend in all three 
directions. Clearly, the only way to fully represent such 
a deformation is to use a fully quadratic element. Any 


element that is less than fully quadratic will be too stiff. 


ПИ LEMPERATURE GRADIENT CASE STUDIES 

The temperature gradient case studies used the same 
models as the flexure case study. Cases 1, 2, 3 refer to 
the same orientation of the 12-node brick as in the flexure 
study. Case 4 is the 20-node brick model. 

Hot pressed silicon nitride was used as the specimen 
material. The coefficient of thermal expansion is given in 
ES [6]. A 100°F temperature gradient was imposed across 
the short transverse (0.125 inch). The stress-free refer- 
ence temperature was 1500%F. A variably-noded (12-21) 
three-dimensional brick from the ADINA finite element pro- 
gram was used to model the problem. Nodal temperatures were 
input via a temperature tape. This procedure is described 


in Appendix A. 
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As in the flexure case study, the normalized deflections 
of the y-z plane of symmetry were plotted versus distance 
along the span. Beam theory was used as the basis for com- 


Epson. The deflection predicted by beam theory zs [7]: 


= Q 2 
τει) κα 
EN. ὁ, Ξ -al‘(T, = T,)/2t 
` А 2 ај: 
@ = eoefficient of thermal expansion (inch/inch/ °F) 
t = depth of beam 
(Т„ - Т.) = temperature difference through depth of beam 


Цена ве валсе between the lower reactions 


Two different case studies evolved from the temperature 
problem. Reference [9] states that the Bernoulli-Euler as- 
Sumptions can be used in practical analyses of beams under 
thermal loadings. The assumptions state, in part, that the 
effects of lateral contraction can be ignored, i.e., Poisson's 
ratio (v) may be taken as equal to zero. In the case studies 
that follow, Poisson's ratio was used in the first study, and 
set equal to zero in the second. 

1. Results of the σα ο, HENO Ss Sons OA) 

The normalized deflection curves for this case study 
are shown in Fig. 12. The displacements from x=0 to x=L are 
plotted. Cases 2 and 3 are stiffer than beam theory. The 
curves show that the maximum deflection calculated for these 
two cases is 75% of the maximum deflection predicted by beam 


theory. 
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The results of Case 4, which used the fully quadratic 
20-node brick, show excellent agreement with simple theory. 

Figure 12 shows that the results of Case 1 (the 12- 
node brick with quadratic side along the span) predict a de- 
flection greater than theory. It is puzzling that the re- 
sults show that the element is even more flexible than the 
20-node brick. The temperature algorithm used by ADINA was 
checked extensively. No errors were detected. No satis- 
factory explanation for the 12-node brick behavior has been 
found. It is suspected that the inherent anisotropy of the 
12-node brick may be the cause of the problem. However, 
extensive further examination would be necessary to resolve 
the dilemma. 

Possibly more interesting than the displacement re- 
Sults are the stresses calculated. Recall that if a linear 
temperature gradient were imposed on a bar that had no ex- 
ternal constraints, the bar would deform into an arc of a 
circle. Also there would be no stress in the bar according 
to a fundamental theorem of elasticity. In this case study, 
the bottom nodes at x = 0.75 inch were fixed in the z-direc- 
tion. The part of the specimen to the left of this con- 
straint (1.e., toward the free end) should, however, show no 
stress until one reaches the immediate vicinity of the con- 
straint. Each of the cases that use the 12-node brick pre- 


dict stresses, both normal and shear, on the order of ΠΠ, 


¡O πο’ PSI in the free end. The 20-node brick, however, 
predicts stresses of approximately 5 PSI or less in this 


region. OE ¡A OS in/in/°F). 
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The above discussion concerning the 12-node brick 
models indicates the possible uncertainty that can result 
from using this element in thermal calculations. 

2. Results of Case Study with v = 0 

This case study was undertaken in an attempt to ex- 
plain the results of the Case 1 orientation of the 12-node 
brick discussed above. Everything was the same as in the 
previous study, except that the effects of Poisson's ratio 
were neglected. The results are shown in Fig. 13. The Case 
2 and 3 models have become significantly stiffer. The maxi- 
mum deflection was approximately 67% of the value predicted 
by theory, as compared to approximately 75% when Poisson's 
ratio was used. The results for Case 4 have remained the 
same (to three significant digits). The displacements for 
Case 1 have decreased to such a degree that they were almost 
identical to Case 4. The spurious free end stresses cited 
previously remained. 

In summary, by ignoring the effects of Poisson's 
ratio, three of the models have become significantly stiffer, 
while the fully quadratic model has remained substantially 
the same and in close agreement with the theory. 

It is interesting to note that if the specimen were 
truly unconstrained and had a linear temperature gradient, 
the results of analyses that take into account Poisson's 
ratio should agree exactly with the results of analyses that 


ПО [9]. 
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А И ЗАО ОНИ ИО О ENG 


This chapter discusses the techniques developed to con- 
struct a finite element mesh for the airfoil section of a 
ceramic gas turbine blade. This particular engine airfoil 
lent itself to detailed analysis because of the relative 
ease with which the mesh could be generated. This is due 
to the linear fairing between the tip and root sections of 
the airfoil. Although linear fairing was not optimal from 
an aerodynamic point of view, it was desirable from a manu- 
facturing point of view. Its relative ease of fabrication 


caused it to be used in practice. 


DE VELOPMENT OF AIRFOIL PROFILES 

The Garrett Corporation provided Mylar cross-sectional 
profiles of the first stage rotor blade of their experimen- 
tal ceramic turbine [10]. These profiles were for various 
radial distances as measured from the axis of rotation 
along the stacking axis of the blade. In addition to an 
outline of the cross-section and its relation to the stack- 
ing axis, a table of x-y coordinates was given for the pro- 
file. These coordinates described 62 points on the profile. 

The 62 profile points were assembled into x and y ma- 
trices. These matrices were used to develop cross-sectional 
plots on the time sharing terminal and later on the Calcomp 


plotter. Figures 14 and 15 show representative cross- 
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sections of the tip and hub sections. The рае points 
are marked by an 'x' and are connected by straight line 
Sie gments. 

Since there was linear fairing between the tip and hub, 
it could be easily verified that the coordinates listed on 
the Mylar drawings were given at the same relative positions 
along the foil profile. In other words, the ith point on 
mae R=3.) inch profile corresponded to the ith point on the 
R=4.0 inch profile, where R is the radial distance along the 
stacking axis from the axis of rotation. 

The root profile was then divided into elements. Care 
was taken to develop elements with convex corners. Also, a 
relatively fine mesh was desirable in the leading and trail- 
ing edge areas of the airfoil to enable accurate modeling 
of the anticipated stress concentrations there. The result- 
ing leading and trailing edge elements proved quite trouble- 
some; this will be discussed more fully below. 

The mesh developed in the airfoil root section dictated 
the mesh in the tip section. This was due to the unique 
correspondence of points in the tip and hub sections. The 
tip section was then connected to the hub section by linear 
nodal coordinate generation.  Connectivities were generated 
from the leading to the trailing edge. This scheme took 
maximum advantage of the automatic generation features of 


SAP IV and ADINA. 


B. THE NEGATIVE JACOBIAN DETERMINANT 
As mentioned earlier, the leading and trailing edge ele- 


ments were ill-behaved. The leading and trailing edges of 
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this blade are arcs of a circle. When an element with qua- 
dratic sides 15 used, these elements are prone to re-entrant 
angles, and resultant negative or zero determinants of the 
Jacobian matrix. In other words, the inverse of the Jacobi- 
an matrix does not exist because there is not a unique cor- 
respondence between the natural and local coordinates of the 
element [8]. This will stop the program because the inverse 
of the Jacobian matrix is used in the strain-displacement 


Matrix [B]. Recall: 


ΤΕ. 
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[81^ (C1 [B] dCVOL) 


Also: [K] - Ivor 
EE - inverse of the Jacobian matrix 
[K] = Stiffness matrix 
{u} = Н] {о. ) - displacement field 
[H] - matrix of shape functions 
[C] So elas tie properties 
le) Spend rax 


r,S,t = natural coordinates 
As the above equations show, the inverse of the Jacobian 
is an essential part of the calculations for the strain ma- 
trix, and the stiffness matrix. Experience has also shown 
that when the determinant of the Jacobian is small Com or 


smaller) the displacements at that node may lose almost all 


their numerical significance. 


MA MESE VERIFICATION THROUGH PSAP1 

Kibler [11] implemented a finite element preprocessor 
for mesh verification for the SAP IV and ADINA programs. 
MaS enables the user to check visually the coordinates and 
geometry of a mesh by generating graphical displays of the 
completed mesh. Oblique orthographic projections are pro- 
duced using the NPS Calcomp model 765 plotter. Other dis- 
plays are possible. 

The PSAP1 program is especially useful because of the 
various mesh presentation options in the program. Among 
them are: 

1. Node numbering 


2. Element numbering 
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ои <p leoded plots 

н. Full rotation about each axis 

5. Sectioned plots based on geometry, 
node number, or element number 

6. Displacement postprocessing using either 
the deformed structure, or displacement 
vectors at the nodes. (See Chapter VII) 

Other options, as well as a User's Manual, are fully dis- 
cussed in Ref. [11]. 

Figures 16 through 19 show the various mesh presentations, 
for the linear airfoil sections.  PSAPl was used extensively 
in the production of the airfoil mesh. The clear graphical 
presentation assures the engineer that he is using a mesh 
free of geometric or connectivity errors. 

The elements with negative Jacobians were separately 
selected from the mesh and greatly magnified. The re-entrant 
angles were clearly visible on these plots as shown in Fig. 
20. Large scale plots of the hub and tip cross-sections for 
these Ec were made by hand. The mid-size nodes were 
moved inward toward the geometrical center of these elements. 

Two airfoil meshes were developed. One was composed of 
20-node bricks, and the other of 12-node bricks. The 12-node 
bricks were oriented such that the quadratic side could be 
used to model the airfoil camber. This orientation was the 
same as the one used by Garrett. Since the 12-node brick 
has only one quadratic side, part of the leading and trailing 


edge curves were eliminated. This is shown in Fig. 21. Aside 


30 





from this difference, the geometry represented by the two 


meshes was identical. The 12-node airfoil is shown in Fig. 


22. 


Пре пете се о: each of these meshes are shown 


en Table I. 
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K ΟΕ ΡΕ STUDY 


After thoroughly investigating the ceramic bar problen, 
the next step was to analyze linear airfoil section of the 
ceramic gas turbine blade using the ADINA code 12-node and 
20-node bricks. The bar problem showed that significant 
differences occurred as a result of the brick type and 
orientation used. The purpose of the airfoil study was to 
determine the difference, if any, between the 12-node brick 
airfoil and a 20-node brick airfoil. Each airfoil was 
loaded using the actual design point conditions that could 
be expected in the ceramic gas turbine. Two general types 


of loading were investigated: centrifugal and thermal. 


A. CENTRIFUGAL LOAD 

A centrifugal load preprocessor was developed at NPS, 
and 1s discussed fully in Ref. [12]. The preprocessor can 
De used to generate consistent nodal loads for a given 
angular speed. The rotor design speed of 41,730 RPM was 
used. 

Figures 23, 24 show the deformed foil superposed on the 
Umndeformed foil. (Chapter VII discusses using PSAP1 to plot 
the deformed structure.) The root plane coincides exactly 
in all cases because the nodes there were fully constrained. 
Two general deformations can be seen: first, there is a 


general elongation of the blade in the radial direction, and 
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second, there is an untwisting of the blade. Figures 25-26 
show the deformed 12-node and 20-node foils superposed on 
each other. It can be seen that the 12-node model is stif- 
fer in representing the untwisting deformation. Tables II 
and III show comparisons of stresses and maximum principal 
stresses, as calculated by the 12-node and 20-node bricks, 
at selected nodes in the root section. At each of the tabu- 
lated nodes, the stresses predicted by the 20-node brick are 
higher than those predicted by the 12-node brick. The 
stress contours for the root section are shown in Figs. 29 


mud 30. 


bee TAERMAL LOAD 

Reference [13] shows the steady state isotherms. For 
this study, the isotherms were modeled as constant for 
various radial distances along the foil. The cited reference 
shows that this is a good assumption in the root section, 
where the isotherms are flat and closely spaced. In the 
upper portion of the foil, however, the isotherms have a 
definite curvature in the radial direction. Figure 31 shows 
the nodal temperatures used in this study. Preparation of 
the nodal temperature tape is discussed fully in Appendix A. 

Figure 32 shows the deformed foil superimposed on the 
original geometry. A general enlargement of the foil in all 
directions is seen. Table IV shows a comparison of thermal 
messes at selected points. (It should be noted that these 
are not evaluated at nodal points, but rather at the Gauss 
points. ADINA does not give nodal stresses for thermal load- 
ss.) 
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eee CONCLUSIONS FROM AIRFOIL STUDY 

This comparative study shows how the choice of a finite 
element can affect the results. The 12-node brick, even 
though oriented in the most advantageous way, 15 stiffer 
than the fully quadratic three-dimensional element. The 
stiffer element produces lower stresses in the root section. 
It is these root section stresses which transmit the tensile 
forces and moments to the blade platform and eventually to 


the critical dovetail area. 
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πο oto le ΕΙ PREPROCESSOR 


The version of ADINA installed at NPS does not contain 
a distributed load feature. Finite element texts, such as 
Ref. [14], often state a consistent allocation of a uniform 
surface load for a rectangular element. The consistent 
load for a curved surface, such as the airfoil studied in 
this report, 1S more complex. Since a distributed load is 
used to model the aerodynamic pressure on an airfoil, a 
program to calculate consistent pressure loads for curved 


surfaces was necessary. 


E DEVELOPMENT OF THEORY 
A consistent load is defined as a set of nodal forces 
which produce the same amount of work as the original force. 


Stated mathematically: 
= 5 
<v.>{u,} = Í ARPA ou c um S AREA) 


where: EE ΛΕΕΙ πο а cons stent Loads 


nodal displacements 


na 
= 
He 
ES, 
II 


«f?» = vector of surface forces 
i Ld IBS - displacement field 


[N;] = shape functions 
Consider the general surface shown in Fig. 33. 


Work = of 


A 


u 
n pdA-u = <v,> Γ᾽ E) 
W 
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n = unit normal to the surface (positive outward). 
p * pressure on dA. 
u = displacement vector = ul + vj + wk 


For ease of computation, it is desirable to perform the 
integration in the natural coordinate system. The vector 


ndA is transformed to the natural coordinate system by: 


98 dR 
ndA = ae x 58) da ав 
where: on) eee E POSi tron vector of point on 


surface as shown in Fig. 33. 
аа ав = differential area in natural coordinate 
system. 
a,8 are dummy variables standing for the appropri- 
ate natural coordinates which define the 


pressure face (See Table V). 


IR 

b This is a vector tangent to a line a where B 
is constant. The components of this vector 
are elements of the Jacobian matrix. 

aR 

55 is similarly defined. 


The pressure and displacement fields can be discretized 


with the same shape functions: 


P = <N;>1p;) (2) 
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а 
V = 0 
W Ü 
where: М. = 
|. 
{р:) = 
ш 
1 = 
W: 
E 
Let 
IR 98 
= X 58? do 
then: 
ndA*u = <A, 
Az 
A, 
Az 


Bsing Баз (2 


1 


1. 

М 0 V ( 
1 1 
0 IM W; 


shape functions 


matrix of nodal pressures 


matrix of nodal displacements 


dB = (Aji + Aoi + Agk) да ав (4) 


u 
A, λα ν да ав 

w 
E τ; να 
HR oov 3,8 
ELE 25a ^1.8 


can now be written 


1 


N; 0 0 
<A AD Ag? 0 N; 0 V <N,>{p,}da ав 
0 : 
-1 3E + 


S 





l fi N. 0 o | 
= | 
E - «ληλρλα7 0 N. O i <N.>{p.] da dg (6) 
АС 
e - 0 0 x 


БИН PROGRAM STRUCTURE 
The consistent pressure load program is structured as 
follows: 
Subroutine TRANS 
- reads element number, coordinates, and connectivity 
- reads the face to which pressure applied (NFACE) 
- reads pressure at each node of NFACE 
- reads number of Gauss points to be used 
in surface integration (2-6) 
= calls CUBAT 
Subroutine CUBAT 
- establishes vector of weights and sampling 
points to be used in Gaussian integration 
- calls SHAPE 
Subroutine SHARE 
- evaluates shape functions at Gauss points 
- evaluates derivative of shape function at 
Gauss point 
- forms Jacobian matrix 
- evaluates determinant of Jacobian matrix to check 
for zero or negative Jacobian determinant 


= сас ЕАСЕР 
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Subroutine FACEP 
- establishes matrices in integrand 
- performs matrix multiplication to establish 
integrand 
- performs surface integration 
A listing of the consistent pressure load program is in 


apendix C. 


С. PROGRAM VERIFICATION 

Several test cases were run to insure the accuracy of 
the pressure load preprocessor. A uniform distributed sur- 
face load was rotated through all faces of a one inch cube, 
and the 2-6 Gauss point integration scheme was checked for 
each face. The results agreed exactly, both in magnitude 
and sign, with the consistent load allocation found in 
her. [114]. 

The program was successfully checked to insure that dis- 
tributed surface force could be applied to more than one 
face of the element simultaneously, and that pressures could 
be applied simultaneously to faces that share a common edge. 

The ceramic flexure specimen mesh (Case 4) was also 
used to verify an actual problem. Pressure loads were 
applied to various faces of the specimen, including adjoin- 
ing faces. In each case, the program produced correct re- 
Sults. 

As a final test, a distributed load was also placed 


across the span between the supports. A sketch of the 
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problem is shown in Fig. 34. 


pared with simple beam theory. 


shown in Fig. 35. 


Pee USER'S MANUAL 


The results were again com- 


The results were excellent 


Figure 36 shows a hexahedral element in the natural co- 


ordinate system. 


оси сп Ејс. 35. 


tion used in the pressure load preprocessor. 


ШОШО corner nodes and four midside nodes. 


tribution acting on an element 
the nodal pressure intensities 
nodes. For example, if Face 4 


Р Р 


00? Pis> Pig» and Pay must be 


The node numbering convention is also 


Table V shows the face numbering conven- 


Each face has 

The pressure dis- 
face is defined by specifying 
at the corner and mid-side 


were loaded P P P | 


ОРИ 


specified in that order. 
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Appendix В lists the input data and results for the sample 


Fu problem, shown in Fig. 36. 
iS Control Card (1015) 

Columns Variable Meaning 

1-5 NUMNP Total number of node 
points in problem. 

6-10 NUMEL Total number of 
elements in problem. 

11-15 NGP Number of Gauss points 
to use in surface 
integration. (2<NGP<6) 

16-20 NCUR Number of ADINA load 


curve that Output loads 


refer to. (usually, NCUR=1) 
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21-25 


26-30 


31-35 


36-40 


1-55 


2. 


11-50 


Nodal Pressures 


NPLOAD 


ΠΠ τ 


TERENDE? 


υπ Κιτ» 


IPRINT 4 


να 1215. 


NEREL (J) 


NFACE (Ј) 


Eo >) 
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Number of element 
faces with distributed 
loads. 

l-program will just 
read and print ADINA 
put deck 

l-program will not 
print element-by- 
element solution. 
l-program will not 
ие ошро deck orf 
consistent loads. 
l-program will omit 


printing input data. 


ΕΕ) 


element number 

face to which pressure 
applied. (See Table V) 
nodal pressure 
intensities input here: 
corner nodes first, 
then midside. (See 
Table V) Positive pres- 
sure is directed in 
same sense as outward 


Uait normal: 





КИНЕ СОЕ PRESSURE LOAD PREPROCESSOR ON CP/CMS 
The pressure load preprocessor is ideally .suited for use 
on CP/CMS. A compiled version of the program (filetype=text) 
should be stored on a private disk or workspace. A data 
file, such as FTO04+F002, should also be defined. The input 
data deck (Appendix B) can then be read into this file [15]. 
It is suggested that an executive type routine also be 
used: 
PED EGER EROR 
τ ΤΡ: υπ GENE JENIOS FODBDPERCPERM) 
БЕРЕБИ Об РЕ (PERM) 
ΕΕ ΕΡΕΕ ο PUN (PERM) 
ο ΕΝ РИМОВ ОБЕ 
© асе а 
This defines your input/output files and executes the 
program. Output 15 directed to the line printer, and card 
punch. 
The pressure program utilizes approximately 300K BYTES 


Of core. The extra core should be requested at LOGIN by: 


Emp SK 
where: 
nnnn = user number 
tt = terminal number 
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И ПА ESA AS CA DISPLACEMENT POSTPROCESSOR 


PSAP1 is capable of postprocessing displacements. The 
results can be presented in two ways: 

1. The deformed structure can be plotted. In this case, 
scaled displacements are added to the original geometry to 
produce the deformed coordinates. The deformed structure is 
then plotted using this new geometry. 

2. The displacements can be represented as scaled vec- 
tors. The scaled displacement vectors are plotted at origi- 
nal geometry nodes. 

The purpose of this chapter is to explain the modifica- 
tion of the ADINA code to store the displacements, and how 


to use the PSAP1 postprocessing capabilities. 


A. MODIFICATION OF ADINA 

Reference [16] states that a version of SAP IV had been 
modified to punch a deck of displacement cards from the out- 
put data. When dealing with large quantities of data, how- 
ever, storage of this data on direct access storage devices 
is preferable. Therefore, the ADINA code was modified to 
Store the output displacements on device 58, in addition to 
Ae usual printed output. 

The flag for saving the displacements is the integer 
number 1 in column 80 of the ADINA master control card [4]. 
The flag variable name is ITP 58, and is passed in common 


ARST. 
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The job control language (JCL) for ADINA must also be 


meditied to define device 58: 


Ep ИО ПОАРО UNTT=3330., VOL=SER= 
[Т DSN=Snnnn.pppp 
where: 


пл = tne user number 


пој 


ΡΡΡΡ is the name assigned by the user to identify 


the data set. 


Sean on ro rmatted Werte Statement is used to store 


the data on device 58, a variable blocksize must be used 


when allocating disk space. Reference [17] 


procedure, 


БИН MODIFICATIONS OF PSAPL 


explains this 


There are presently two subroutines in PSAP1 that can 


be used to input displacement data to PSAPl. 
1. Subroutine DATA 9 was added to PSAP1 
This subroutine reads displacement data from 
Reference [16] discusses its use, and should 
2. Subroutine DATA 5 reads displacement 
from device 58. A JCL card must be added to 


pxdefine device 58: 


Бу bos] 
punched cards. 
be consulted. 
data directiy 


ea l] 


//G0.FT58001 DD DISP=0LD,UNIT=3330,VOL=SER=DISKO, 


m DSN-Snnnn.pppp 


USER'S MANUAL FOR DISPLACEMENT POSTPROCESSING 


The PSAP1 User's Manual [11] applies with the following 


meet) cations: 


a) Namelist Option 
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Variable, Value 


NUDISP 


NVDISP 


NWDISP 


KDATA 


NVALUS 


1 


al 


nnn 


py Namelist Pict 


KDISP = 


ТОМАС 


DMAGS = 


n 


τ 


Description 
x direction displacements input 
y direction displacements input 
z direction displacements input 
subroutine DATA 5 used to read 
displacement data from device 58. 
nnn is the integer number of dis- 
Јеша ments set, to) be plotted. 
This is equal to the number of 


node points. 


n=l: plot of deformed structure 
n=3: displacements represented 
by vectors at ... 
n=1: direct magnification of 
displacement data by DMAGS. 
n=2: scaling of displacement 
data to a maximum value 
of DMAGS. 
ris a floating point number 
which indicates the magnifica- 


tion of the displacements. (0.1 


has been found appropriate.) 


Figures 37 and 38 are examples of displacement postpro- 


cessing. 


The structure is the bending specimen loaded in 


four-point flexure. 
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Figure 37 illustrates the deformed structure option of 
PSAP1. Figure 38 shows the displacements of element 15 


plotted as scaled vectors at the node points. 
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VIII. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

1. The 12-node brick is a stiffer model of a structure 
than the fully quadratic 20-node brick. 

2. The accuracy of the results obtained using a 12-node 
brick depends on geometry, loading, and the orientation of 
mae element. 

3; The flexure problem showed that two orientations of 
the 12-node brick are stiffer than the third. Resulting dis- 
placements and stresses differed significantly, and in all 
cases were lower than the 20-node brick and theory. 

4, The airfoil study showed that the maximum principal 
stresses predicted in the root section could differ by as 
much as 14.9% depending on the basic element used. The 12- 
node brick under-predicted the 20-node results. The use of 


the 12-node brick could therefore produce an unsafe design. 


В. RECOMMENDATIONS FOR FUTURE WORK 

As is true with many projects, this thesis raised many 
new questions. The following areas of future work are 
Exccested: 

l. The ADINA code should be reviewed with an eye toward 
Streamlining it. A basic code using only two- and three- 
dimensional elements, and the linear isotropic material 


model would be a useful tool. A reduced version of ADINA 
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covld then be installed on the time-sharing system, or 
other stand-alone system. 

2. The Field program, developed at NPS, should be ex- 
amined and expanded if necessary. This program would serve 
as an excellent thermal preprocessor. It should be modi- 
fied to produce a nodal temperature tape that could be 
used in ADINA, or other code capable of using temperature 
Ырп}. 

3. Aerodynamic pressure loads should be examined. А 
program that could calculate the nodal pressure intensi- 
ties for an element would provide the necessary input in- 
formation for the distributed load preprocessor. 

4. Centrifugal, pressure and thermal loadings should 


ре superposed. 
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Figure 3. Two-Dimensional Plane Element 
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SCHEMATIC SHOWING SELF-ALIGNING CAPABILITY 


Figure 4. Schematic of 4-Point 
Flexure Rig 


92 





τ E 


IL 0.75" 1.125" 
/ 
0.75" 1.5" | 0,75" 
R R 
А B 


Figure 5a. Bend Specimen Free-Body Diagram 


р 


Figure 5b. Bend Specimen Shear Diagram 
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Figure 5c. Bend Specimen Moment Diagram 
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Figure 6b. x-z View of Half-Span with Guided 
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Figure 6c. y-z View of Specimen with Guided 
Midplane Boundary Condition 
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PLANE OF 
ME TRY 


Figure 7. Cross-Sectional View of Deformed 
Specimen Showing x~z Plane of 
Symmetry 
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Figure 8a. 12-Node Brick Orientation for Case-l 





Figure 8b. 12-Node Brick Orientation for Case-2 





Figure 8c. 12-Node Brick Orientation for Case-3 





Figure 8d. 20-Node Brick; Case-4 
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Figure 9b. Free Body Diagram for Temperature Gradient 
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Airfoil Mesh Exploded Plot with 


Elements Numbered 


Figure 19, 








Figure 20. Trailing Edge Element Showing 
Re- entrant Angle 
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Figure 21. 12-Node Airfoil Section Showing 
Truncation of Leading Edge Element 
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ΕΞ: CHARACTERISTICS 


20-NODE BRICKS 
ου DE  _... ss coe ccc cco sce es ces eee LU=node brick 
Pe OL NOdES..crcccccccscccccvevesevetos 
Нел ог ејетепт а... с... с... с с с с. с «ес 99 2 
Шштсеос5 ОГ Ггегдот. с . 2 . с... 2... .........1152 
Пи = | Бапаы:асһ.....................93 
Кит nalt=-bandwidth.....oo.oooooooo..129 


Mimber of stiffness matrix elementSs.....107064 


IRZNODERBRIECKS 
ПОИНТА Ое .........................12-подае brick 
ПИ појегз......... с... «се. с«се«•««• 270 
ИДО ејепепт 5... ...... he o. eee. со со сео «02 
Mes OL тгесдош. . . 2... .. со. «<< « « .... e DO 
ELM  г—рапаилаЕћ,. о с с с с о ШБ 
Цитати half-bandwidth........ eee eee 66 


Number of stiffness matrix elements.....29106 
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TABLE II. COMPARISON OF MAXIMUM STRESSES IN AIRFOIL ROOT 
SCIELO (Ог, Ете. 27) 


O xx yy 022 
LOCATION POIL CEST) CRS IE, (PST) 


I 1 2553 2553 10213 
20 2658 2658 10634 
B ilt? 2897 2897 11588 
20 ρα Sa 12843 
C 1 3082 3082 oo 
20 3560 3560 14245 
D 12 3130 3130 12520 
20 3656 3656 14626 
E | i? Sg 3074 12300 
20 3625 3625 15000 
Ε 12 2946 2946 11784 
20 3447 3447 13788 
8 12 2781 2781 {ШЧ ЛЕ 
20 893 ΞΟ η 12924 
H 152 2608 2608 10432 
20 3000 3000 12001 
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Ш ВЕЕ III. 





COMPARISON OF MAXIMUM PRINCIPAL STRESSES IN 
ΠΗ, ος Εεσο (cf. Fig. 27) 


Ae 


97, (20) - 9, (12) 


0 Difference = crema 
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ο n ТУ. COMPARISON OF SELECTED THERMAL STRESSES IN 
ΝΕΟΙ ΠΠ SECTION (cf. Fig. 28) 
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ПИЛЕ pA PREPARATION OF TEMPERATURE TAPE 


If temperature input is used with the ADINA program, 
nodal point temperatures are stored in sequence on device 
number 56 [4]. At NPS, ITEL 3330 disks have been used, 
although any storage device may be designated. NSTE + 1 
records are read, where NSTE is the number of solution 
steps. (In linear static analysis, NSTE usually equals 1.) 

The read statement in the program is: 

READ (56) TIME, (TEMPV(I), I=1, NUMNP). 

On each record, the problem time of the nodal tempera- 
tures is read, followed by an implicit "do-loop" from 1 to 
the total number of node points (NUMNP) to read each nodal 
point temperature. 

A program to establish the temperature tape is listed 
at the end of this Appendix. It will be discussed in three 
parts: (A) job control language, (B) the RDADIN subroutine 


and (C) the binary write. 


A. TEMPERATURE TAPE JOB CONTROL LANGUAGE (JCL) 

The program listing in Appendix B starts with a two step 
procedure [17]. The first is a purge operation that returns 
the previously reserved disk space to the operating system, 
the second step reserves disk space for the data set which 
is named F0099.TAPE. One cylinder of 3330 disk space (248K 


bytes) is reserved on one of the NPS computer center resident 


° 





disks. The data set has a variable block size, and an ex- 
piration date of 360 days after establishment. 
The JCL employed combines ease of continued use of the 


disk space, with minimum chances of data set expiration. 


ΤΠΕ RDADIN SUBROUTINE 

The RDADIN subroutine was originally programmed as part 
of the PSAP1l preprocessor. The subroutine will read an 
ADINA input deck from title card through element connectivity. 
It generates node numbers, nodal coordinates, element numbers 
and element connectivities for the entire mesh. Each of the 
above is stored in a linear array. The programmer is. then 
able to utilize any or all of these arrays to input the 
nodal point temperature. 

For example, in the linear airfoil mesh there were iso- 
thermal planes. The matrix of z-coordinates (radial distance 
from axis of rotation) was used. Each z-coordinate was 
checked by a series of 'if statements', and nodal tempera- 


tures were assigned. 


Pome lif, BINARY WRITE 

After the matrix of nodal temperatures was established, 
the next step was to write it sequentially onto device 56. 
Array RCD1 contained the problem time for the temperatures, 
and temperatures for each node; i.e., it was a linear array 
of dimension NUMNP * 1. This array could then be written 
onto device 56 by the statement: 


ИЕ Coo) RCD I. 
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Device 56 must then be defined in the GO step by: 


εὐ. ΕΠ5ΡΕΏΠΙ DD UNIT=3330,VOL=SER-DISKO1l, 
fy ~WON=FO099.TAPE,DISP=OLD 


Each time the WRITE statement is executed one record is 
established. The sequence must be repeated until NSTE * 1 
records are established, where NSTE is the number of solu- 
tion steps. 

In the linear airfoil problem, there was only one solu- 
tion time step. Therefore, only two records on device 56 
were needed. One record was at time=0, and the next record 
was One time step increment later. 

It is recommended that a hard copy of the RCD1 matrix 
be printed to ascertain the accuracy of the temperature 


каре. 


55 





20099 . ТАРЕ, РОКБСЕ 


DISKO1,DISP=SHR 


JOB CARD) 


A 
VOL=SER 


3330=DISKO1 ,DSNAME 


Or ο) 
шс — 
XA = 
WNN AX 
ar et De 
мс Әс 


APPENDIX B. 


SER=DISKO1 »DISP=(NEWeKEEP), 


330,VOL= 
1,1,1)) 


3 
( 


? 


4 
L 


> 


> 


ш 
xa 
шс 

LL. 


NNNSNN ж МММ ММ. 
~ 


“МММ 


0-2) 
3K >K 2 >K ας XK Ἂς “κ Ὡς Ἂς ας ἃς XK 3K 3⁄2 SK XK X K 3 X EEK FH KK FR FT KERN 


+ 


end CO UI) L PUSO 
ао £ <t << Z Ú2 
TNT |Ο 
e О ЈУ OD 
ο ΖΞ D 
Ш о сИ. ОМ 
CQ KZ eo 
с) О we DW 
= GOE ° > 
αι) ΣΣ Го 
O Ш с“ -e 
wu i Z O m U pe 
О Оши 10 
Die „л e Z 
wnuxn ~ 
UJ eil] mls ect 
О Ш202 002 > 2 
Or ДР о — <£ 
Zog- € 
JO =Z LL > >< 


2 ed. Xa 

Ce OA 
T OQ «c «ou 
> ст гл 


55 


NODAL TEMPERATURES ARE WRITTEN SEQUENTIALLY CN DEVICE 58 
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NODAL TEMPERATURE TAPE PROGRAM 
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